library(tidyverse)
library(sf)
library(terra)
library(openxlsx)
library(mapview)
library(RColorBrewer)
library(osrm)
library(movecost)
library(fs)
library(stringr)
library(DT)
library(synthdid)
library(summarytools)
mapview::mapviewOptions(fgb=FALSE)
load("data/china.RData")
load("data/sakhalin2.RData")
load("data/sakhalin3.RData")

衛星夜間光データとは

ntl1992<-terra::rast(x="data/9828827/Harmonized_DN_NTL_1992_calDMSP.tif")
ntl2018<-terra::rast(x="data/9828827/Harmonized_DN_NTL_2018_simVIIRS.tif")
ntl1992 <- ntl1992 %>%
  terra::crop(y=ext(c(123.1,146.9,24.0,45.7)))
ntl2018 <- ntl2018 %>%
  terra::crop(y=ext(c(123.1,146.9,24.0,45.7)))

1992年日本

terra::plot(ntl1992)

2018年光り輝く日本　若干赤くなってる？

terra::plot(ntl2018)

### What is NTL(night-time lights) ・夜間光は近年民生・産業部門での複合的な経済指標の代理変数として社会科学領域で活用が広まっています

・NTL(夜間光)データについては、中小発展途上国のさらに細かい行政区において経済変数としての意味を成すかという議論が活発に行われている。

・Are night-time lights a good proxy of economic activity in rural areas in middle and low-income countries? Examining the empirical evidence from Colombiaは、コロンビアを例に夜間光と詳細な経済指数を比較し分析した結果、都市部の方が精度は高いが、農村部でも正の相関が確認されており、非常に有用なデータであることを示している。

・特に、経済状況を示すデータにアクセスできないような地域（今回示す北方領土やチベット地域）を時系列で分析する際に有用と考えることができる。

・今回用いるのは、1992-2018年までのNTLであり、Harmonized VIIRS & DMSPという二種類の計測方法を重み図けしてつなぎ合わせたものを利用する。

Synthetic DID法とは

発想

因果推論アプローチの一つであるPotential outcome frameworkにおいて重要なのは 観察しえない反実仮想\(Y_{itD=0}\)をいかにして導き出すかという問題です。

反実仮想というのは例えば、統一した東ドイツ（に当たる地域の）経済\(Y_{itD=1}\)に対する統一しなかった場合の東ドイツ経済\(Y_{itD=０}\)や拓銀が倒産した後の北海道経済（観察可能）に対する拓銀が生存した場合の北海道経済（観察不可能）などです。

今回紹介するSDID法は、本当に簡単に言うと

拓銀が倒産しなかった世界線の北海道を、他のいろんな県を組み合わせて作ろう！

という発想が元になっています。

説明

SDID法は差分の差分法(DID法)と合成コントロール法(SC)を組み合わせた手法です。

DID推定量は以下の通りです （書き方は様々あるが）\[ (\hat{\tau}^{DID}, \hat{\mu}, \hat{\alpha}, \hat{\beta})={\arg\min}_{\mu,\alpha,\beta,\tau} \sum_{i=1}^N \sum_{t=1}^T(Y_{it}-\mu-\alpha_i-\beta_t-D_{it}\tau)^2 \] 各ユニット・時間における結果変数\(Y_{it}\)は定数\(\mu\)、時間効果\(\beta\)、個人効果\(\alpha\)の線形結合で表されます。今回求めたい因果効果をとらえるのは処置変数\(D_{it}\)の係数\(\hat{\tau}^{DID}\)です。このモデルの最小化を行うことでDID推定量は求められます。

一方、合成コントロール法はSynthetic Control Methods for Comparative Case Studies(Abadie et al. 2010) で提案されました。

アイデアとしては、処置を受ける前まで処置群を統制群（処置を受けない群）の重み付け（加重平均）で近似できるような重み \(\hat\omega\) を求めて、その加重平均で処置後の処置群を予測し、その差分を処置効果とするものです。それを表す式は以下の通りです。

\[ \hat{\tau}_{it}^{synth}=Y_{Nt}-\sum_{i=1}^{N-1}\omega_{i}Y_{it} \]

問題の重み\(\hat\omega_{i}\)は処置前までの期間、つまり\(T-1\)期に注目して以下のL2正則化項（罰則付き）二乗誤差最小化問題で求められます。

\[ \hat{\omega}_{i}=\underset{\omega}{\arg\min}\frac{1}{T-1} \sum_{t=1}^{T-1}(\sum_{i=1}^{N-1} \omega_{i} Y_{i t}-Y_{N t})^{2}+ \frac{1}{2}\zeta\|\omega\|_2 \]

…つまり、どちらも処置群と統制群のバイアス（元の違い）と処置前と処置後のバイアス（放っておいても変化するかも）を統制して差分(処置効果）を取りたい訳です。そのなかで、いかに”並行トレンド仮定”という非常につよい仮定を緩められるかという所が問題になっていますね。

ですが、合成コントロール法ユニットごとの重み\(\hat\omega_{i}\)しか使っていません。もちろん、重要な年と重要じゃない年を峻別するような時間の重みも考えられますよね!

そのような時間に対する重み\(\lambda_{t}\)を合わせてDID推定量の中に組み込んでしまったのが今回使用するSynthDID推定量です。(\(\hat\lambda_{t}\)の推定に関しては\(\hat\omega_{i}\)の最小化問題の対象がユニットから時間に変わっただけです。)

\[ (\hat{\tau}^{SDID}, \hat{\mu}, \hat{\alpha}, \hat{\beta})={\arg\min}_{\mu,\alpha,\beta,\tau} (\sum_{i=1}^N \sum_{t=1}^T(Y_{it}-\mu-\alpha_i-\beta_t-D_{it}\tau)^2 \hat{\omega}^{SDID}_i \hat{\lambda}^{SDID}_t) \]

うれしさ

・簡単に扱える

・説明・解釈がしやすい

・少ない処置群＆大量の対照群

ロシア「クリル(北方領土)開発計画」効果検証

「クリル諸島(北方領土)社会経済発展プログラム (2007-2015)」とは…

ーそれだけ、ロシア政府は東部前線基地の開発計画に投資しているのだ。

このプログラムには180億ルーブルの資金が投入され、前任者とは異なり、全額が融資されます。 新しいインフラの構築と既存のインフラの近代化に重点が置かれます。
交通機関- 岸壁、空港、道路
エネルギー - 送電網と新しい発電能力
社会的なもの - 病院や幼稚園の建設、新しい住宅の建設など…

-連邦政府目標プログラム「クリル諸島の社会的・経済的開発」より

mapview(RUS %>% filter(ID_1==60))

クリル諸島(北方領土)は1996年にロシアによる開発計画が策定されましたが、国内経済等の混乱もあり立ち消えになっていましたが、2007年から本格的な資本投下が行われ、主に民生部門での発展が計画されました。 しかし当該地域は特殊な状況に置かれており、どの程度発展したのか分析できるようなデータにアクセスすることは出来ません。

しかし夜の光を隠すことは出来ないでしょう。

Descriptive statistics

データ前処理（省略）

ntl<-function(x){
  i<-terra::extract(
    x=x,
    y,
    fun="sum",
    na.rm=T
  )
  i<-select(i,!ID)
  return(i)
}

RUS<-sf::read_sf(dsn="data/RUS_adm/RUS_adm3.shp")

y<-terra::vect(sf::read_sf(dsn="data/RUS_adm/RUS_adm3.shp"))

#mapview(RUS %>% filter(NAME_1 %in% c("Kamchatka", "Khabarovsk" ,"Sakhalin")))

russia <- lapply(filecon, ntl)#ほんとに時間かかります

rus1<-bind_cols(russia)

rus1<-rus1 |> 
  as.tibble() |> 
  rename_with(\(x) str_replace(x, "Harmonized_DN_NTL_", " "))

rusdata<-bind_cols(RUS,rus1)

rusdata<-rusdata |> 
  as.data.frame() |> 
  select(!c("ID_0","ISO","NAME_0","TYPE_3","ENGTYPE_3","VARNAME_3","geometry","NAME_3","NL_NAME_3"))

rusdf<-rusdata |> 
  pivot_longer(cols = c(` 1992_calDMSP`: ` 2018_simVIIRS`),
               names_to = "year",
               values_to = "ntl") |> 
  mutate(lntl=log(ntl+1))

rusarea<-expanse(y, unit="m", transform=TRUE)/1000000

rusarea<-rusarea |> 
  as.tibble() |> 
  mutate(ID=1:length(rusarea))

rusdf2<-left_join(rusdf,rusarea,by=c("ID_3"="ID"),copy=T)

rusdf2<-rusdf2 %>% 
  mutate(perarea=ntl/value,
         lperarea=log(perarea+1))

rusdf2$year<-as.double(str_sub(rusdf2$year,2,5))

east<-rusdf2 %>% 
  filter(ID_1 %in% c(3,9,83,82,24,56,28,60,40,61,   12))

極東連邦管区（北方領土も含む）の夜間光(対数)/面積データ

使用したデータは後述します

DT::datatable(east %>% select(NAME_1,NAME_2,year,lperarea))
summary(east %>% select(ntl,lntl,value,perarea,lperarea))
##       ntl                lntl            value             perarea         
##  Min.   :       0   Min.   : 0.000   Min.   :     0.2   Min.   :  0.00000  
##  1st Qu.:    1011   1st Qu.: 6.919   1st Qu.:  3033.7   1st Qu.:  0.07882  
##  Median :    3680   Median : 8.211   Median :  9825.0   Median :  0.45290  
##  Mean   :   29452   Mean   : 7.636   Mean   : 29476.7   Mean   :  7.05588  
##  3rd Qu.:    8456   3rd Qu.: 9.043   3rd Qu.: 34714.7   3rd Qu.:  2.11444  
##  Max.   :11078183   Max.   :16.220   Max.   :320110.2   Max.   :295.14318  
##     lperarea      
##  Min.   :0.00000  
##  1st Qu.:0.07587  
##  Median :0.37356  
##  Mean   :0.88117  
##  3rd Qu.:1.13605  
##  Max.   :5.69084
east<-rusdf2 %>% 
  filter(ID_1 %in% c(3,9,83,82,24,56,28,60,40,61,   12))

クリル諸島(北方領土)の夜間光の伸び(対数）

kril %>% 
  ggplot() + 
  geom_line(aes(x = year, y = lperarea))

NTL analysis

east<-east %>%  #データ整形
  filter(year<=2015) %>% 
  mutate(treatment=if_else(ID_3==1820,if_else(year>=2007,1,0),0)) %>% 
  as.data.frame()

setup = panel.matrices(east,time="year",treatment= "treatment",
                       outcome ="lperarea",unit="ID_3") #計画行列

tau.hat = synthdid_estimate(setup$Y, setup$N0, setup$T0) #推定

top.controls = synthdid_controls(tau.hat)[1:10, , drop=FALSE] #ウェイト計算（上位10個)

se = sqrt(vcov(tau.hat, method='placebo')) #std.erro

点推定・95%信頼区間

sprintf('point estimate: %1.2f', tau.hat)
## [1] "point estimate: -0.01"
sprintf('95%% CI (%1.2f, %1.2f)', tau.hat - 1.96 * se, tau.hat + 1.96 * se)
## [1] "95% CI (-0.46, 0.45)"

推定結果(プロット)

synthdid::synthdid_plot(tau.hat)

青の折れ線はクリル諸島（北方領土）のデータであり、赤のグラフは統制群の加重平均によって近似された反実仮想である。

synthdid_units_plot(tau.hat, units = rownames(top.controls))

重み\(\hat\omega\)の大きさプロット

結果

北方領土の夜間光は仮に社会経済発展プログラムを実行しなかった場合に比べ、平均的に-0.01(log y)つまり1%減少し、効果がある政策とは言えないという結果になった。

しかし他のロシア極東地域も同様に著しく発展していると理解することも可能であり、連邦政府がサハリン２を含む 極東開発に非常に力を入れているのが確認出来る。

今後、中国のゼロコロナ政策撤廃によりエネルギー需要がさらにひっ迫することが予見されており、西側諸国に制裁を受けているロシアは極東のエネルギー開発に活路を見出す可能性がある。

中国青海-西蔵鉄道開通による経済効果分析

工事中です。

データ・参考文献

  1. Synthetic Difference In Differences(Arkhangelsky et. al., 2019)を読んだhttps://dropout009.hatenablog.com/entry/2019/05/06/185620

  2. その無茶振り，（Rで）GISが解決しますhttps://rpubs.com/k_takano/r_de_gis